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We investigate the (conservative) dynamics of binary black holes using the Hamiltonian formula- 
tion of the post-Newtonian (PN) equations of motion. The Hamiltonian we use includes spin-orbit 
coupling, spin-spin coupling, and mass monopole/spin- induced quadrupole interaction terms. We 
investigate the qualitative effects of these terms on the orbits; in the case of both quasicircular and 
eccentric orbits, we search for the presence of chaos (using the method of Lyapunov exponents) 
for a large variety of initial conditions. For quasicircular orbits, we find no chaotic behavior for 
black holes with total mass 10-40 Mq when initially at a separation corresponding to a Newtonian 
gravitational- wave (GW) frequency less than ~ 150 Hz. Only for rather small initial radial distances 
(corresponding to a GW frequency larger than ~ 150 Hz), for which spin-spin induced oscillations 
in the radial separation are rather important, do we find chaotic solutions, and even then they 
are rare. Moreover, these chaotic quasicircular orbits are of questionable astrophysical significance, 
since they originate from direct parametrization of the equations of motion rather than from widely 
separated binaries evolving to small separations under gravitational radiation reaction. In the case 
of highly eccentric orbits, which for ground-based interferometers are not astrophysically favored, 
we again find chaotic solutions, but only at pericenters so small that higher order PN corrections, 
especially higher spin PN corrections, should also be taken into account. Taken together, our surveys 
of quasicircular and eccentric orbits find chaos only for orbits that are either of dubious astrophys- 
ical interest for ground-based interferometers or which violate the approximations required for the 
equations of motion to be physically valid at the post-Newtonian order considered. 
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I. INTRODUCTION 

Relativistic binary systems made of compact objects, 
such as neutron stars or black holes, are among the most 
promising (or, at least, are among the best-understood) 
candidates for the production of gravitational waves de- 
tectable by both ground- and space-based gravitational 
wave observatories. The difficulty of detecting the signals 
from such systems has led to a theoretical effort to un- 
derstand the gravitational waveforms likely to be emitted 
by such systems, which in turn has led to a considera- 
tion of their dynamical behavior. In particular, several 
authors Q, |^ y, ^ |^ have investigated the presence of 
chaos in the dynamics of compact binaries, motivated in 
part by the effect of such chaos on the calculation of theo- 
retical templates for use in matched filters. The extreme 
sensitivity on initial conditions that characterizes chaotic 
systems would lead to significant difficulties in the imple- 
mentation of such filters, since the number of filters would 
grow exponentially with increasing detection sensitivity. 

In the extreme mass-ratio limit, chaos was found in the 
equations describing a spinning particle orbiting a non- 
rotating (Schwarzschild) black hole Refs. 0,01 ex- 
tended this result to the case of a rotating (Kerr) black 
hole, finding widespread chaotic solutions. As demon- 
strated in 0, however, the values of the total spin for 
the test particle leading to chaotic solutions are not real- 



izable in physical systems. Furthermore, |y| showed that 
chaos, while widespread for these unrealistic spin values, 
disappears in all cases for physically realistic spins. In 
short, there is strong evidence that extreme mass-ratio 
systems, which are most relevant for proposed space- 
based gravitational wave detectors, are not chaotic for 
any parameter values of physical interest. 

The case of the comparable-mass binaries more rel- 
evant to ground-based gravitational-wave observatories 
has been investigated by several authors H, Q using 
the post-Newtonian (PN) equations of motion in the La- 
grangian formalism, using harmonic gauge 0, • There 
was initially some doubt regarding the results presented 
in 2], which found chaos in the PN equations for spin- 
ning bodies, since the timescale of the chaos was not re- 
ported: it was not clear whether the chaos discovered in 
the equations — in the conservative limit neglecting gravi- 
tational radiation reaction — would have time to manifest 
itself in the inspiral timescale tinsp- Furthermore, the 
work in Q cast doubt on the presence of chaos in these 
systems, finding that the Lyapunov characteristic expo- 
nents for the PN equations, which measure the divergence 
rate of nearby trajectories, are zero in all cases tested. 
However, Q found some initial conditions, correspond- 
ing to rather eccentric orbits, that do have positive Lya- 
punov exponents, indicating the presence of chaos, with 
characteristic times shorter than tjnsp, raising the possi- 
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bility that theoretical templates calculated for systems 
with spinning compact objects are affected by chaos. 

In the present study, we examine and extend these re- 
sults by investigating the dynamics of spinning binary 
black holes using a Hamiltonian formulation of the post- 
Newtonian equations of motion ITfll Ull in the 
Arnowitt-Deser-Misner (ADM) gauge. In order to make 
chaos formally possible, we exclude gravitational radia- 
tion reaction; since tests for chaos technically require an 
infinite-time limit, the finite inspiral times due to radia- 
tion reaction would eliminate the possibility of chaos. On 
the other hand, we do include post-Newtonian terms in- 
volving the spin of the two bodies: the addition of spin is 
essential to create the possibility of chaos, since without 
spin the constants of the motion constrain the motion to 
be at most quasiperiodic. As discussed in Sec. Ill Al we 
use four separate spin terms in the equations of motion to 
model accurately their effect of the dynamics. Wc focus 
on black holes, to the exclusion of other compact astro- 
physical objects, because two of these spin terms (which 
involve spin quadrupole effects) are known exactly only 
for black holes, and yet their magnitudes are comparable 
to the spin-spin coupling and hence cannot be ignored. 
(This is an extension of previous work, as other authors 
have not considered these quadrupole terms when inves- 
tigating chaos.) 

In Sec.ni'ws write down the PN Hamiltonian, including 
spin terms, and in Sec. lIIII we discuss how we choose ini- 
tial conditions for quasicircular and eccentric orbits. We 
then investigate chaos for comparable-mass binary black 
holes (Sec. IIV|I . Since binary black hole inspirals tend 
to circularize under gravitational radiation reaction, we 
focus first on the important special case of quasicircular 
orbits and then analyze eccentric orbits. As in previous 
work, we favor Lyapunov exponents fSec. lIV'X|) to quan- 
tify the presence (or absence) of chaos. 

We work almost exclusively in geometric units (G = 
0=1). Euclidean vectors, such as appear in the post- 
Newtonian equations of motion, are set in boldface, and 
we use vector arrows to denote relativistic 4- vectors. The 
symbol log refers in all cases to the natural logarithm. 

II. THE POST-NEWTONIAN EQUATIONS OF 
MOTION 

The post-Newtonian (PN) equations for the two-body 
problem are an approximation to full general relativ- 
ity, essentially involving a series expansion in v/c. As 
in the case of the classical two-body problem, in the 
post-Newtonian case it is possible to describe the mo- 
tion of a relativistic binary in the ccntcr-of-mass frame. 
A typical orbit is shown in Fig. ^ In this paper, 
we use the center-of-mass Hamiltonian, as developed in 
Refs. H m m El m. The Hamiltonian formulation 
is particularly convenient for our present purposes: since 
detecting chaos involves determining the separation of 
nearby phase-space trajectories, it is convenient to work 



directly in terms of spatial coordinates and their corre- 
sponding conjugate momenta — a criterion automatically 
satisfied by the Hamiltonian formulation. 

A. The Hamiltonian formulation 

We can represent the PN Hamiltonian schematically 
as follows: 

H = Hjsi + iJpN + Hso + Hss- (2.1) 

We include the following terms: Newtonian, post- 
Newtonian (usually through 2PN order, i.e., v*/c'^, and 
sometimes through 3PN order, i.e., v^/c^), spin-orbit 
coupling (through 1.5PN order), and spin-spin coupling 
(through 2PN order). (We omit the radiation reaction, 
as discussed in the introduction.) Throughout this treat- 
ment, we use X for the (relative) position, P for the con- 
jugate (relative) momentum, and (81,82) for the spins 
of the two objects. 

We denote with mi and m2 the mass of the two bodies, 
and introduce the total mass M = nii + m2 and the re- 
duced mass /i = mim2/M} Using these mass variables, 
we can express the first term (the standard Newtonian 
energy) as follows: 

H.^^-^. (2.2) 

2^ r 

where r = |X|. The post-Newtonian terms used in this 
paper are IPN and 2PN (and, in some sections, 3PN), 
given by i E [H E IH : 

HpN = /i(i/iPN + H2PN + -^spn), (2.3) 

where 

i^iPN = i(3r7-l)(p2)2 

^hs + rj)p' + rjin-pf]- + ^ (2.4) 
2 q 2(7^ 

and 

16 

+ i[(5-2077-3r;2)(p2)2 

-2,7'(n-p)V-3,7^n.p)4]l 

Q 

+ i[(5-K8ry)(p2)-H3,7(n.p)2]l 
z q 

-i(l + 3r,)l (2.5) 

and 
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(a) (b) 

FIG. 1: The orbit of two maximally spinning 10 Mq black holes, using a Hamiltonian formulation of the post-Newtonian 
equations of motion, (a) The orbit embedded in Euclidean space; (b) the projection onto the x-y plane. Lengths are measured 
in terms of the total mass M — mi + m2, and we show a schematic horizon at rn = M, indicating the collapse radius where 
the relative separation of the two bodies is the sum of their horizon radii. (We have rn,i = rn^, and rii,2 = m2 for maximally 
spinning black holes, so the relative separation of collapse is rn = mi + m2 = M.) Note that, in contrast to Newtonian orbits, 
the orbit is not closed, and the orbital plane precesses around the center of mass. 



For clarity in what follows, we adopt the arbitrary convention 
that mi > m2. 



^3PN (q, P) = (-5 + 357? - IQif + iW) (P')' 

+ 1 [(-7 + 42,7 - 5377^ - Sr;^) (p2)3 + (2 _ 3r;)r;2(n • p)2(p2)2 + 3(1 _ n)r,\n ■ p) V - ■ vf] - 

16 q 



^ (-27 + 136ry + 10^) {p^f + 1(17 + iQiM^ ■ P)'p' + + 43ry)ry(n ■ p)^ 
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(2.6) 



In the preceding formulas, we have r/ = mim2/M'^, and 
we use the unit vector n = X/r and the reduced canon- 
ical variables p = P//i and q = X/M. In most of this 
paper, we will measure momenta in terms of fi and dis- 
tances in terms of M, so we will typically not distinguish 
between the canonical and reduced canonical variables. 
In this paper, in TJpN we use only the terms through 
the 2PN corrections (thereby omitting the 3PN terms) 
except where explicitly noted. 

The next term in Eq. 1)2.1(1 is the spin-orbit coupling 
(corresponding to Lcnsc-Thirring precession in the ex- 
treme mass-ratio limit m2 ^ mi): 



L-S. 



so 



(2.7) 



where 

Se.= f2+|l^)si+f2 + |l^)s.. (2.8) 
\ 2 mi J \ 2 7712 / 

The spin coupling term in Eq. 1(2. 1(1 has three compo- 
nents: 

H-ss = i?SiS2 + -f^SiSi + Hs^Si- (2-9) 
The first term, the spin-spin coupling, is 

Hs^s, = ^[3(Si • n)(S2 • n) - Si • S2], (2.10) 

which is valid for all bodies (e.g., neutron stars or white 
dwarfs). The next two terms we include are monopole- 
quadrupole interaction terms, and their form is valid only 
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for black holes.^ They are [11113 



HsA = 7^ [3(Si • n)(Si • n) - Si • Si] ^ (2.11) 

2r'^ mi 



and 



1 777 1 

^s.S. = ^ [3(S2 . n)(S2 ■ n) - S2 • S2] ^. (2.12) 

With the full Hamiltonian in hand, we can now derive 
the equations of motion using Poisson brackets. As in 
classical Hamiltonian mechanics, the time-evolution of a 
dynamical quantity /(X, P, Si, S2) is simply the Poisson 
bracket of the quantity with the Hamiltonian: 



(2.13) 



The Hamiltonian equations of motion for the (relative) 
position and (relative) momentum are then the familiar 
canonical equations: 

dx_ m 

To derive the spin equations of motion, we use the canon- 
ical angular momentum Poisson bracket 



which yields 



and 



dSi _ dH 



dS2 _ dH 



X Si = rji X Si 



X S2 = r22 X S2. 



(2.15) 



(2.16) 



(2.17) 



Eqs. H2.14|l and (|2.16|I - H2.17|I . with the Hamiltonian given 
by Eq. (|2.1|l , are the equations of motion used throughout 
this paper. 



B. Conserved quantities 

There are many conserved quantities in the post- 
Newtonian equations. These constants of the motion con- 
strain the dynamical behavior of the system and provide 
valuable checks when testing a numerical implementation 
of the equations. Here we discuss all the quantities known 
to be conserved, and at which orders they are conserved. 



1. Quantities conserved at all orders 

The following quantities are conserved at all orders: 

• Total energy H: H = {H, H} = by the antisym- 
metry of Poisson brackets. 

• Total angular momentum J = L + S1 + S2: see |15| . 

• The spin magnitudes ^i and 52: this is evident 
from Eqs. (|2.16|) and (|2.17|) . since the cross product 
is perpendicular to the spin and hence can change 
only its direction. 

2. Quantities conserved only through spin-orbit coupling 

If we neglect the terms quadratic in the spin (i.e., we 
include only terms through spin-orbit coupling), the fol- 
lowing additional quantities are conserved: 

• = L ■ L: at this order, L — Scs x L/r^, which 
changes the direction of L but not its magnitude. 



• L • Sf, 



see '15!. 



In our numerical implementation, we verify that the 
above quantities are conserved at the proper orders — that 
is, we check that energy, total angular momentum, and 
the spin magnitudes are always conserved, and verify that 
L and L • Scff are conserved when including only terms 
through spin-orbit coupling. To calculate the Lyapunov 
exponents we use the techniques and routines described 
m Chapter 5 of [3. For more details on our numerical 
implementation, see the Appendix. 



III. PARAMETERIZING POST-NEWTONIAN 
ORBITS 

We discuss here two convenient methods for parame- 
terizing initial conditions for post-Newtonian orbits. We 
first describe a method that gives orbits that approx- 
imately satisfy desired values of eccentricity, pericenter, 
and orbital inclination. We then treat the important spe- 
cial case of quasicircular orbits. Finally, we examine the 
effects of varying the post-Newtonian terms included in 
the Hamiltonian. 



A. Eccentric orbits 

One convenient method for fixing initial conditions 
starts with the eccentricity e and pericenter r^, which 
allows for the easy creation of bound orbits and makes 
contact with the Newtonian limit. For a Newtonian 
orbit with given values of (e,rp), we can calculate the 
corresponding (reduced) angular momentum Ln and en- 
ergy Ejq using 



^ We leave the generalization to neutron stars and other compact 
bodies to future work. 



Ln — 



(3.1) 
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(a) (b) 

FIG. 2: A highly eccentric orbit for two maximally spinning 10 Mq black holes, (a) The orbit embedded in Euclidean space; 
(b) the projection onto the x-y plane. Lengths are measured in terms of the total mass M = mi +m2, and we show a schematic 
horizon at rn = M. The energy and angular momentum of the orbit correspond to a Newtonian orbit with eccentricity e = 0.9, 
pericenter rp = 3.7 M. The empirical values of these parameters, as determined from the numerical solution to the equations 
of motion, are e = 0.899 and rp = 3.7 M. The spins are Si = 0, 1/^2) and S2 = (1/2, 0, V3/2). 



and 

= 1 - (3.2) 
2rp 

where we have included the rest energy in the energy 
term. Our parameterization method then involves find- 
ing a post-Newtonian orbit with the same energy E as 
the Newtonian orbit, with the canonical momentum 
set equal to L^.^ 

The details of the eccentricity parameterization involve 
first setting Xq = {rp, 0, 0) and pr = 0. We then set p^ = 
LTss/r (the Newtonian value of the (j) momentum), leaving 
only pg undetermined. At this point we could set pg = 0, 
in agreement with the Newtonian value; this choice leads 
to perfectly valid initial conditions. But, since we are 
studying the dynamics of a Hamiltonian system, we wish 
to assign a privileged role to the energy, and the choice 
pg — does not lead to a post-Newtonian system with 
the energy calculated from Eq. H3.2|l . Hence, we force the 
energies to agree using 

H -^Newtonian 5 

(3.3) 

which gives a sixth-degree polynomial equation in pg. 
Since the Newtonian value of pg is exactly zero, in or- 
der to produce the post-Newtonian orbit analogous to 
its Newtonian counterpart we choose the real root of 
Eq. 1)3. 3|l closest to zero. (This choice still leaves two 
roots, corresponding to initial values of pg in the ±z di- 
rection. We arbitrarily choose the negative root, so pg is 
initially in the +z direction.) 



We use lower-case for (the canonical momentum conjugate to 
</>) to distinguish it from = ■ P, the component of P. 



We should note that the conditions Xq = {rp, 0,0) and 
Pr — 0, which are chosen for computational convenience, 
mean that the initial orbit has its pericenter in the x-y 
plane, but this condition is not true for a generic post- 
Newtonian orbit; as a result, the initial conditions do 
not necessarily satisfy the requirements for a valid post- 
Newtonian orbit (particularly for the nonrotating case 
in which orbits are confined to a plane). Fortunately, 
as soon as the orbit reaches the true pericenter — that 
is, when pr is again — those initial conditions do re- 
sult in a valid post-Newtonian orbit. (This means that 
the radius of the pericenter requested — i.e., rp — and the 
empirical pericenter differ, but by examining the numer- 
ical solutions we find that the requested and empirical 
pericenters typically differ by only a few percent.) Since 
no finite segment of the orbit can affect the presence of 
chaos, which is defined as an asymptotic property of the 
system (Sec. IIVX|) . the small invalid piece at the begin- 
ning of the orbit does not affect the final result. 

As a result of this parameterization method, we are 
able to find solutions to the post-Newtonian equations 
of motion that have empirical values of e and rp (as de- 
termined by examining the numerical solution directly) 
quite close to the values of corresponding Newtonian or- 
bits (Fig. 01. 



B. Quasicircular orbits 

A second, more specialized parameterization of the FN 
initial conditions enforces the condition of quasicircular- 
ity. In particular, through third post-Newtonian order 
the quasicircular orbits are in fact exactly circular, and 
even with spin-orbit coupling added there exist "spheri- 
cal" orbits, i.e., orbits confined to lie on a sphere, with 
fixed radius but varying angle 9. Once any of the spin- 
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spin terms is turned on, exact sphericity is impossible 
in general, but it is still possible to satisfy exactly the 
conditions leading to spherical orbits in the absence of 
spin-spin coupling. These orbits are especially important 
for modeling possible sources of gravitational radiation, 
since the orbits of compact binaries are expected to cir- 
cularize due to gravitational radiation reaction 

The conditions leading to quasicircular orbits are as 
follows. Given an initial radius rg, we set (/)o = and 
9q = 7r/2, so that 

Xo = (ro,0,0). (3.4) 
We then require that the initial radial momentum vanish: 



orbits to make the effects especially easy to see, but this 
is not a necessary restriction. Figs. ^ and [3 show that the 
PN equations capture the essential aspects of relativistic 
orbits, such as the characteristic precession of the orbital 
plane. 

Figs. I^HS] show quasicircular orbits satisfying 
Eqs. (|3.4|| - ||3.7|l . with increasing orders (up to 2PN) of 
the PN Hamiltonian included. The third figure includes 
all the spin-spin terms (S'iS'2, SiSi, and S2S2); we do 
not show any orbits with subsets of these terms activated 
because, to the precision visible in the figure, the S'iS'2 
term dominates, and the figure is identical with either 
of the other terms removed. 



{Pr)0 = 0. 



(3.5) 



[Since the Hamiltonian is quadratic in Pr, this means 
that (r)o = {dH/dPr)p^=o — 0, so that (at least initially) 
the radius is not changing.] Finally, we require that the 
initial values of Pr and vanish, which means (using 
Hamilton's equations) that 



and 



dPr 
~dt 



dH 
dr 



[dPe 



= 



0. 



(3.6) 



(3.7) 



Given the initial position and the initial spins, these 
equations can be solved numerically for the initial val- 
ues of Pg and P^, thereby giving a complete set of initial 
conditions. 

When the spin-spin terms arc included, we can always 
initially satisfy the conditions for quasicircularity [Eqs. 
(|3.5|) ~ (|3.7|l ] , but these conditions are not preserved by the 
evolution, since in this case orbits with constant radial 
separation no longer exist; the radial position oscillates as 
a function of time. These spin-induced radial oscillations 
can have non-negligible amplitudes at small separations, 
and it is unclear whether they represent orbits that could 
result from the adiabatic inspiral of quasicircular orbits 
under radiation reaction. A better method would be to 
set quasicircular initial conditions when the black holes 
are rather far apart (so that spin effects are negligible) 
and then evolve the system toward the final plunge by 
including radiation reaction in the dynamics. The cal- 
culation of radiation reaction effects for the Hamiltonian 
framework with spin couplings is currently under com- 
pletion |T3|; we plan to include them in the dynamics 
and investigate chaotic behaviors in the near future. 



C. Post-Newtonian orbits at various orders 

Here we show some of the effects of turning on or off 
the various post-Newtonian terms. We use quasicircular 



IV. INVESTIGATING CHAOS IN THE 
POST-NEWTONIAN EQUATIONS 

Previous studies of chaos in the post-Newtonian equa- 
tions considered coniparable mass-ratio binaries with ec- 
centric orbits 0, 0, 113, IS 112 ■ Here, we also consider 
comparable mass-ratio binaries and focus first on quasi- 
circular orbits and then on eccentric orbits. As noted in 
the introduction, these orbits are particularly important 
because many astrophysically relevant binary systems 
should circularize due to the energy lost to gravitational 
radiation. We set the radii of these orbits so that the 
frequencies of their corresponding gravitational waves lie 
in a range 40 Hz < /gw < 240 Hz, roughly correspond- 
ing to the frequency band for the Laser Interferometer 
Gravitational- Wave Observatory (LIGO) and VIRGO. 
[More properly, we choose radii such that the Newto- 
nian frequency of the gravitational waves (which is sim- 
ply twice the orbital frequency) lies in the LIGO /VIRGO 
band, as discussed in Sec. lIVBl below.] For initial separa- 
tions such that the Newtonian GW frequency is smaller 
than 40 Hz, we never find chaotic quasicircular orbits. 



A. A brief discussion of Lyapunov exponents 

As in previous works Q , Lyapunov exponents are 
our primary tool for investigating the nonlinear dynam- 
ics of general relativistic systems. We have discussed at 
length in our techniques for calculating these ex- 

ponents for systems similar to the PN equations. Here 
we present a brief summary of Lyapunov exponents. 

Given an initial condition in the phase space of a dy- 
namical system with n degrees of freedom, we imagine an 
n-dimensional ball of nearby initial conditions centered 
on that point. As the dynamics unfold, in general the ball 
is stretched in some directions and squeezed in others, de- 
forming into an n-dimensional ellipsoid under the action 
of the flow. Such an ellipsoid has n principal (semi) axes, 
and the average rate of stretching or squeezing of each 
axis is a Lyapunov number, whose natural logarithm is 
the Lyapunov exponent associated with the axis. In gen- 
eral, the zth Lyapunov exponent of a dynamical system 
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FIG. 3: A post-Newtonian quasicircular orbit of two maximally spinning 10 M© black holes, w ith Newtonian, IPN, and 2PN 
terms turned on. (a) The orbit embedded in Euclidean space; (b) the radius r = ^/aJ^ + jp + z^ as a function of time. Lengths 
are measured in terms of the total mass M — mi + m2, and we show a schematic horizon at rn = M. 




FIG. 4: A post-Newtonian quasicircular orbit of two maximally spinning 1 Mq black h oles, with terms through spin-orbit 
coupling, (a) The orbit embedded in Euclidean space; (b) the radius r = ^yx^+y'^ + z^ as a function of time. Lengths are 
measured in terms of the total mass M = mi -I- m2, and we show a schematic horizon at ru = M. The addition of spin-orbit 
coupling to the N, IPN, and 2PN terms destroys exact circularity, but exact sphericity is preserved. 



is 

A. = li„ M, ,4.1, 

t^oc t 

where ri{t) is the ith principal ellipsoid axis. Implement- 
ing this prescription numerically leads to a visualization 
of the exponents as a plot of log [ri{t)] vs. t (so that the 
slope is the exponent Ai), which we refer to as a Lya- 
punov plot. The result for a chaotic PN orbit appears in 
Figs. El and El 

In practice, following the evolution of the phase-space 
ellipsoid, and thereby extracting all the Lyapunov expo- 
nents of the system, involves using the Jacobian matrix 
of the system to model an "infinitesimar' ball that cap- 
tures the true linear approximation to the dynamics. It 
is also possible, and computationally faster, to extract 
only the largest Lyapunov exponent by considering only 
one nearby initial condition, joined by some small devi- 
ation vector to the original point. In what follows, most 
of our simulations use this faster (but less robust) devi- 
ation vector method, but we have checked many of the 
results using the Jacobian method. Further details of the 
various techniques for calculating Lyapunov exponents 
appear in Q and especially [l^. 



It is worth noting that we can think of the PN system 
as constrained^ since we wish to think of the spin mag- 
nitudes as fixed. In other words, given an initial spin 
vector, a "nearby" initial spin should point in a differ- 
ent direction but have the same magnitude. The system 
thus has only ten true degrees of freedom (three for rela- 
tive position and momentum, and two for each spin) , and 
should therefore have only ten Lyapunov exponents. The 
constraints lead to significant complications in calculat- 
ing the Lyapunov exponents; see (16] for several methods 
of addressing these complications. 

The principal value of the largest Lyapunov exponent 
is that it provides the e- folding timescale t\ = l/\ for the 
divergence of nearby trajectories. The formal definition 
of A in Eq. (|4.1|l requires an infinite-time limit, but of 
course any numerical method for A must introduce some 
finite cutoff. As a result, in general it is impossible to 
say with any certainty that a system is not chaotic — 
even if it appears that Amax for some ^cutoff, chaos 
may yet manifest itself on longer timescales. Neverthe- 
less, it is possible to calculate nonchaotic baseline orbits 
(corresponding, for example, to the PN terms through 
the spin-orbit coupling), whose nearby initial conditions 
still exhibit some nonzero (power-law) separation. If a 
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FIG. 5: A post-Newtonian quasicircular orbit of two maximally spinning 10 Mq black holes, with all the terms from Sec 111 Al 
(including the all the spin-spin couplings) present, (a) The orbit embedded in Euclidean space; (b) the radius r = ^^/aP^+ip^+z^ 
as a function of time. Lengths are measured in terms of the total mass M = mi -I- m2, and we show a schematic horizon at 
rn = M . The addition of spin-spin coupling to the N, IPN, 2PN, and spin-orbit terms destroys exact sphericity, but Eqs. 13.411 - 
13.711 are still satisfied, leading to nearly circular orbits for these initial conditions. 





2500 5000 7500 10000 12500 15000 
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FIG. 6: The orbit of two maximally spinning 10 Mq black 
holes. The dynamics are chaotic, as shown in Fig. Q The 
initial conditions satisfy the requirements for quasicircular- 
ity, though in fact the orbit's radius is not even approx- 
imately constant (see Sec. II V B 31 below). The initial ra- 
dius is 5.658 M, corresponding to an orbit with a Newtonian 
gravitational- wave frequency of /gw = 240 Hz [Eq. H4.4|l ]. 
The initial spins are Si = (0.13036,0.262852, -0.955989) 
and S2 = (0.118966, -0.13459, -0.983734) ml. The other ini- 
tial conditions are fixed by the conditions for quasicircularity 

fsec. inra . 

suspected chaotic orbit has a Lyapunov exponent with a 
magnitude similar to a baseline orbit, we say that it is 
indistinguishable from or consistent with zero, and hence 
is probably not chaotic. 

For the present problem, the most relevant timescale is 
the inspiral time due to the energy loss from gravitational 
radiation. For the quasicircular orbits considered below, 
this is approximately given by the following formula |23| , 

_ 5 / r \4 

^'"^P " 256 ~ Vm J ' ^^■^> 

where M is the total mass and /i is the reduced mass. 
For eccentric orbits we use Eq. (5.14) of We then 



FIG. 7: The natural logarithms of the ellipsoid axes vs. t 
for the system shown in Fig. |S] The slopes of the lines are the 
Lyapunov exponents. Two nonzero exponents are clearly vis- 
ible (A — ±3.2 X 10~^ M~^), but aU the others are consistent 
with zero. The curve with the largest slope corresponds to 
the upper hue in Fig. |H1 There is an apparent ±A symmetry: 
for each exponent -|-A, there is a corresponding exponent —A; 
even the zero exponents approach zero symmetrically. This 
behavior is a characteristic of Hamiltonian systems 12] . 



adopt the criterion t\ < Unsp as an operational definition 
of chaos, which is equivalent to the condition 

Atinsp > 1 condition for chaotic orbit. (4.3) 

On the other hand, if tx > tinsp, then, even if the system 
is formally chaotic in the conservative limit we consider 
here, the chaos will not have time to manifest itself before 
the final plunge. 



B. A survey of quasicircular orbits 

In this section, we elucidate the effects of varying the 
parameters in the PN equations of motion on the pres- 
ence of chaos in the resulting dynamics. In many of the 
examples, we parametrize the orbits by their radii, or. 
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equivalently, by the "gravitational wave frequency" : 

where we restore the factor of G so that the result is 
in Hz. It is essential to note that Eq. H4.4|l is the New- 
tonian gravitational wave frequency, which is valid only 
for radii that satisfy r 3> M . Nevertheless, Eq. H4.4|l pro- 



vides a convenient way to parameterize the initial con- 
ditions by radius in a way that has transparent physical 
significance in the nonrelativistic limit. When we refer 
below to an orbit with gravitational wave frequency of 
(say) 240 Hz, we mean an orbit with a radius that sat- 
isfies Eq. |jl3I) when /^^* = 240 Hz. It is important to 
remember that this is not in general the true frequency 
of the gravitational wave — for example, at 2PN order, for 
equatorial orbits, averaging over an orbit yields |^ l2^l25| 



/•2PN 

Jgw 



IT \ r-^ J 

31 1 

^4^M4 



1/2 



1 1. o ^GM 



ttL • Spff I — 

2 Af 2 \ r 



(Si-S2)-3(L-Si)(L-S2) 



GA/V 



3/2 



- 3 + -?7 + -77^ 



(4.5) 



with L = L/|L|. For a binary with total mass (10 + 
10) M© at radius r = 5.5658 Af (/^^* = 240 Hz), 
Eq. (lOJ gives /g™ ^ {igi^ 194^ 212} Hz when spins are 
aligned with (orbital) angular momentum, zero, and anti- 
aligned with angular momentum, respectively. For the 
same spin orientations at r = 5.742 M (/qw ' = 150 Hz), 
Eq. g3J) gives = {122, 127, 134} Hz. 

1. A chaotic quasicircular orbit 

We find that post-Newtonian orbits satisfying the con- 
ditions for quasicircularity (Sec. IIII Bp can be chaotic 
(though only for rather small initial radial separations), 
as shown in Figs. EHHl Note from Fig. [7|that the Lya- 
punov exponents come in ±A pairs, a characteristic of 
Hamiltonian dynamical systems. Note also that the 
principal exponent calculated using the deviation vector 
method (Fig. |SJ) agrees closely with the largest exponent 
determined from Fig.[71 which uses the more complicated 
Jacobian method to find the exponents. 

2. Varying spin directions 

We illustrate the effect of varying the spin directions 
by generating a large number of quasicircular orbits with 
randomly oriented (maximal) spins. For each spin con- 
figuration, we choose the radius corresponding to a grav- 
itational wave with frequency /q^v'* ~ ^40 Hz (the high 
end of the LIGO/ VIRGO frequency band). The choice of 
radius is motivated by two main factors. First, choosing 
the lowest possible radius (consistent with the abilities to 
detect the corresponding gravitational waves) likely rep- 
resents a worst-case scenario for chaos, since low-radius 
regions correspond to stronger nonlinearities in the equa- 
tions of motion (as noted in 0). Second, minimizing 
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t(M) 



FIG. 8: The natural logarithm of the principal ellipsoid axis 
vs. time for the system shown in Fig. El using the deviation 
vector method. The slope of the line is the Lyapunov ex- 
ponent, which is approximately A = 3.2 x 10~^ using 
a least-squares fit (which agrees closely with the value from 
Fig. which uses the more sophisticated Jacobian method 
to find the exponent). This corresponds to a Lyapunov (e- 
folding) timescale of = 1/A = 3.1 x 10^ M, which is less 
than a fifth of the inspiral timescale. The simulation data for 
a nonchaotic orbit (light) is shown for reference. 



the radius minimizes the inspiral timescale tjnsp, which 
in turn minimizes the computational cost of a final inte- 
gration time significantly longer than finsp- This allows 
us to achieve a better bound on the suspected zero Lya- 
punov exponent, and increases our confidence that ap- 
parent nonzero exponents represent genuine chaotic be- 
havior. In what follows, the final integration time is ten 
times the inspiral time of each orbit. 

We consider the following mass configurations: (20 -I- 
5)A/0, (10+5)A/q, {5+5)Mq, (lO-HlO)AfQ, (20+20)A'/q, 
(20 + lO)Af0, and (15 + 5)Mq. The result of choosing 
N = 500 randomly oriented initial spins for each case 
appears in Tabled We find the presence of chaotic orbits 
for the (10 + 10)Mq and (20 + 10)Mq cases, but we find 
no chaos for any other configuration. Fig. |S1 shows a 
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TABLE I: The prevalence of chaos in post-Newtonian quasicircular orbits at 240 Hz, for spin directions chosen randomly on 
a unit sphere. We calculate the fraction of orbits whose e- folding times t\ = 1/A are less than the inspiral time tinsp, which is 
our operational definition of chaos. The final integration time is ten times the inspiral time. We also include 95% confidence 
intervals for the reported fractions, and we show the average value of A measured in units of the inverse inspiral time for the 
(20 + 1O)M0 configuration (the only case in our simulation with more than one chaotic orbit). The simulation data represent 
500 randomly chosen initial spin directions for each configuration, with the initial radius fixed by requiring a gravitational wave 
frequency of 240 Hz (as determined by the Newtonian formula). 



Configuration 


Fraction chaotic 


95% confidence interval 


Average chaotic Atinsp 


(20 + 5)M0 





[0, 0.00738] 




(10 + 5)M0 





[0, 0.00738] 




(5 + 5)M0 





[0, 0.00738] 




(10 + 1O)M0 


0.002 


[5.06 X 10"^ 0.0111] 




(20 + 2O)M0 





[0, 0.00738] 




(20 + 1O)M0 


0.104 


[0.0777, 0.136] 


1.45 


(15 + 5)M0 





[0, 0.00738] 
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FIG. 9: Lyapunov exponents for 500 quasicircular orbits as a 
function of total angular momentum J for the (10 + 1O)M0 
configuration. The spin for each body is maximal with ran- 
dom initial spin angles {9,(j)) (which overweights the poles). 
The initial radius corresponds to a gravitational wave fre- 
quency of /gw'' ~ 240 Hz. The Lyapunov exponents are 
measured in terms of the inverse inspiral time 1/finsp, so that 
Afinsp > 1 indicates that nearby trajectories diverge by a fac- 
tor of e on a timescale shorter than the inspiral timescale. 
There are 14 such chaotic initial conditions out of the 500 or- 
bits considered for this configuration. They are clustered at 
the low end of the angular momentum range, indicating that 
the spins are closely aligned with each other and anti-aligned 
with the orbital angular momentum. 
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FIG. 10: Lyapunov exponents for 500 quasicircular orbits as 
a function of total angular momentum J for the (20 -|- 1O)M0 
configuration. The spin for each body is maximal and ran- 
domly oriented, and the initial radius corresponds to a grav- 
itational wave frequency of /g^* = 240 Hz. The Lyapunov 
exponents are measured in terms of the inverse inspiral time 
1/tinsp, so that Afinsp > 1 iudicatcs that nearby trajectories 
diverge by a factor of e on a timescale shorter than the inspi- 
ral timescale. There are 49 such chaotic initial conditions out 
of the 500 orbits considered for this configuration. Unlike the 
(10 + 1O)M0 case shown in Fig. |2l the chaotic orbits in this 
case correspond to total angular momentum in the middle of 
the range. 



Lyapunov plot for the strongest chaos in our simulation 
data. The onset of chaos is marked by a transition from 
linear (or at most power-law) separation of nearby initial 
conditions to exponential separation. On our Lyapunov 
plot (which is logarithmic on its vertical axis), chaotic 
orbits appear as linear growth. 

Our initial simulation for the (10 -I- 10)Mq configura- 
tion used random spin angles (0, 0), which does not cor- 
respond to a random orientation but rather overweights 
the poles. This was a stroke of good luck: as shown in 
Fig. El for (1O-|-1O)M0 the chaotic orbits are clustered at 
the lowest values of the total angular momentum J, cor- 
responding to initial spin vectors nearly anti-aligned with 



the orbital angular momentum L (so that J = |L + S| 
is minimized). When running the simulation again using 
randomly oriented spins, we find only one chaotic orbit 
for this configuration. The association of chaos with low 
values of J holds also for the sole chaotic orbit found in 
the (10 -I- 5)Mo case, but it is not a general result: as 
Fig. mi shows, chaos for the (20 -I- 10)Mq configuration 
occurs mainly for values of J in the middle of the possible 
range. 
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FIG. 11: Lyapunov exponents as a function of the (Newto- 
nian) gravitational-wave frequency for quasicircular orbits of 
two maximally spinning 10 Mq black holes. The Lyapunov 
exponents are measured in terms of the inverse inspiral time 
1/iinsp, and the frequencies are chosen to correspond closely 
to the LIGO/VIRGO frequency band. There is an abrupt 
transition to chaos at approximately 160 Hz. 



FIG. 12: The natural logarithms of the principal ellipsoid 
axes vs. time for frequencies on opposite sides of the tran- 
sition to chaos shown in Fig. [TTl /gw" = 158 Hz (light) 
and /g^' = 160 Hz (dark). The slope of each line is the 
Lyapunov exponent, with A ~ 7 x 10~* (light, non- 

chaotic/consistent with zero) and A ~ 4 x 10~^ (dark, 
chaotic). The orbits corresponding to these two frequencies 
appear in Figs. 1131 and 1141 



3. Varying initial frequencies 

We now investigate the results of varying the initial 
(gravitational-wave) frequencies for the strongest chaotic 
orbit from the previous section (as illustrated in Fig. 
The result appears in Fig. ^2 which shows that, even for 
this worst-case scenario, i.e., the system with the largest 
Lyapunov exponent, chaos is absent for initial radii cor- 
responding to /q^' less than 160 Hz. Above 160 Hz, 
there is an abrupt change in the dynamics from regular to 
chaotic (Fig. I12|l , with a maximum Lyapunov exponent 
more than 18 times the inverse inspiral time (meaning 
that nearby trajectories diverge by a factor of e in a time 

tx ~ iinsp/18). 

Cornish and Levin have recently pointed out |26j that if 
unstable orbits are perturbed they could become locus of 
chaos. The abrupt transition from nonchaotic to chaotic 
behavior that we observe in Fig. 1111 when we increase the 
GW frequency (and hence lower the radial separation), 
could correspond to the transition from stable to unsta- 
ble orbits in the spinning ADM Hamiltonian. Thus, to 
better understand the onset of chaos, it would be worth- 
while to explore this intuition further by applying some 
stability criterion in the Hamiltonian framework of the 
kind worked out in the non-spinning case by '^2^ . 

In this particular case, the qualitative change in the 
dynamical behavior from nonchaotic to chaotic is mir- 
rored in the orbits themselves. In particular, the onset 
of chaos is associated with a breakdown in the quasicir- 
cularity of the orbit. As shown in Fig. just below 
the transition to chaos the orbit is nearly circular. Just 
above the transition, despite satisfying the conditions for 
quasicircularity at initial time, the orbits are not even ap- 
proximately circular along the evolution because of the 
strong spin coupling, as shown in Fig. 1141 

To verify that the disappearance of chaos at lower 
frequencies is generic, we repeated the 240 Hz survey 
for the lower end of the frequency range considered 
(/gvk* ~ 40 Hz). The inspiral times are very long in 




FIG. 13: The nonchaotic quasicircular orbit of two maximally 
spinning 10 Mq black holes, corresponding to a gravitational 
wave frequency of /gw'' = 158 Hz. The spins are the same 
as in Fig. |S| The corresponding Lyapunov plot is shown in 
Fig. 1121 The orbit's radius is approximately constant, as re- 
quired for a true quasicircular orbit. 




FIG. 14: The chaotic quasicircular orbit of two maximally 
spinning 10 Mq black holes, corresponding to a gravitational 
wave frequency of /gw'' = 160 Hz. The spins are the same 
as in Fig. 15] The corresponding Lyapunov plot is shown in 
Fig. 1121 Note that the quasicircularity has broken down; the 
radius is not even approximately constant. This qualitative 
change in the orbit accompanies the onset of chaos as the 
frequency increases (with a corresponding decrease in radius), 
as illustrated in Fig. 1111 



12 



TABLE II: The prevalence of chaos in post-Newtonian qua- 
sicircular orbits at 40 Hz, for spin directions chosen randomly 
on a unit sphere. We calculate the fraction of orbits whose 
e-folding times tx = l/\ are less than the inspiral time iinsp, 
which is our operational definition of chaos. The final integra- 
tion time is ten times the inspiral time. We also include 95% 
confidence intervals for the reported fractions. The simulation 
data represent 500 randomly chosen initial spin directions for 
each configuration, with the initial radius fixed by requiring 
a gravitational wave frequency of 40 Hz (as determined by 
the Newtonian formula). In this case, the number of chaotic 
orbits for each configuration is zero. 



Configuration 


Fraction chaotic 


95% confidence interval 


(20 + 5)M0 





[0, 0.00738] 


(10 + 5)M0 





[0, 0.00738] 


(5 + 5)M0 





[0, 0.00738] 


(10 + 1O)M0 





[0, 0.00738] 


(20 + 20)Mq 





[0, 0.00738] 


(20 + 1O)M0 





[0, 0.00738] 


(15 + 5)M0 





[0, 0.00738] 



this case, requiring patience on the part of the simulator, 
but the results are gratifying: as shown in Table HH we 
found not even one orbit with a Lyapunov time less than 
the inspiral time at 40 Hz. Any chaos, if present, mani- 
fests itself in this case on timescales longer than ijnsp- 



4- Varying spin magnitudes 

In Sec. IIV B 31 we created a one-parameter family of 
orbits by taking the worst offender from Sec. IIV 1^21 and 
varying the frequency (or, equivalently, the quasicircular 
radius). In this section, we do the same, but fix the 
frequency and vary the magnitude of one of the spins. 

We showed in Sec. IIV ir2l that the worst-offender orbit 
(with mi = 7712 = 10 M©) is chaotic when 52 = 1 (mea- 
sured in units of TO2). As shown in Fig. 1151 the dynamics 
are nonchaotic for most values of S2, with a transition 
to chaos at approximately S = 0.85. Since the dynamics 
are nonchaotic when S2 = 0, the chaos must be produced 
by the spin terms in the Hamiltonian. 

In Fig. 1161 we show Lyapunov plots for orbits on either 
side of the chaotic transition. Although the difference 
is not as dramatic as the frequency-induced transition 
(Fig. 1121) , there is still a qualitative change in the value 
of the principal Lyapunov exponent. Unlike the system in 
Sec. IIV B 21 this transition does not give rise to a qualita- 
tive change in the orbit as the spin is varied. Instead, the 
chaos manifests itself in the time-evolution of the spins, 
as shown in Figs. 1171 andll8l 
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FIG. 15: Lyapunov exponents as a function of spin for quasi- 
circular orbits of two 10 Mq black holes. We fix the spin of one 
hole at the maximum value (5*1 = ml), and also fix the spin 
directions (which are the same as in Fig. |^ , and then vary 
the spin S2 of the second body. The Lyapunov exponents are 
measured in terms of the inverse inspiral time 1/finsp. There 
is an abrupt transition to chaos when S2 exceeds 0.85 (mea- 
sured in units of m^). 
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FIG. 16: The natural logarithms of the principal ellipsoid axes 
vs. time for values of ^2 on opposite sides of the transition 
to chaos shown in Fig. 115! 5*2 = 0.83 (light) and S2 ~ 0.86 
(dark). The slope of each line is the Lyapunov exponent, 
with A ~ 2.0 X 10~^ (light, nonchaotic/consistent with 

zero) and A ^ 3.5 x 10~^ (dark, chaotic). The orbits 

corresponding to these two spins appear in Figs. I17l and ll8l 



5. Varying the PN terms 

The previous results in this section included all the 
PN terms described in Sec. Ill Al with the exception of 
the 3PN term. Here we investigate the effect of varying 
these terms. We first consider the effect of the 3PN term, 
and then investigate the effects of turning off one or more 
of the spin terms. 

In order to evaluate the effect of the 3PN term on 
chaos in the post-Newtonian equations, we re-calculate 
the simulation of 500 orbits in the (20 + 1O)M0 system 
shown in Fig.^J the result appears in Fig.^1 There are 
some minor differences, including a slight suppression of 
the most chaotic orbit, but the effect of 3PN is not strong. 
This result holds also for the (10 -I- 1O)M0; in this case 
we examine the transition to chaos shown in Fig. ^2 As 
seen in Fig. 1201 the transition occurs at a slightly lower 
frequency, but the difference between the two cases is 
small. Our general conclusion is that the third post- 
Newtonian term has only minor effects on the presence 
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FIG. 18; (a) The chaotic quasicircular orbit of two 10 M© black holes, with spins magnitudes of 5*1 = 1 and S2 = 0.86; (b) 
Cartesian "spin space" showing the time-evolution of Sx, Sy, and Sz. The corresponding Lyapunov plot is shown in Fig. Ilfcil 
There appears to be no qualitative difference between the orbit (a) and the orbit in Fig. llTf a'). but there is a qualitative change 
in the spin behavior. Unlike the frequency transition shown in Fig. 1121 which gives a qualitative change in the evolution of the 
spatial variables, the spin transition to chaos manifests primarily itself in the spin degrees of freedom. 
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FIG. 19: Lyapunov exponents for 500 quasicircular orbits as 
a function of total angular momentum J for the (20 + lO)A/0 
configuration with 3PN term added. The initial conditions 
are identical to those in Fig. 1101 the prevalence of chaos is not 
strongly affected by the presence of the 3PN term. 



of chaos. 

In order to investigate the effect of the spin terms in 
the post-Newtonian Hamiltonian, we focus first on the 
strongest chaos found in our simulations of the (20 -I- 
1O)M0 configuration for /^™* = 40 Hz with 3PN turned 



FIG. 20: Lyapunov exponents as a function of (Newtonian) 
gravitational-wave frequency for quasicircular orbits of two 
maximally spinning 10 A^q black holes with 3PN term added. 
The Lyapunov exponents are measured in terms of the in- 
verse inspiral time l/finsp, and the frequencies are chosen to 
correspond closely to the LIGO frequency band. There is 
an abrupt transition to chaos at approximately 155 Hz. The 
initial conditions are identical to those in Fig. 1111 



on, as shown in Fig. 1101 The Lyapunov plots for a variety 
of PN term combinations appears in Fig. 1211 The most 
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FIG. 21: The natural logarithms of the longest ellipsoid axis r vs. t for the strongest chaos found in the (20 + 10)Mq 
configuration (as illustrated in Fig. 1101 , for varying post- Newtonian terms: (a) 3PN turned off; (b) 3PN turned on. The final 
time is 50 times the inspiral time, and the slopes of the lines are the Lyapunov exponents. In all cases, the Newtonian, IPN, 
2PN, and spin-orbit terms are turned on, but some of the others may be turned off; from top to bottom (with colors, visible in 
electronic versions of this paper, noted parenthetically): all PN terms (black); SiSi and 5'2S'2 turned off (cyan); Si Si turned off 
(blue, nonchaotic); S2S2 turned off (red, nonchaotic); SiSi, S2S2, and S1S2 turned off (orange, nonchaotic). For this particular 
case, the S1S2 spin-coupling term is necessary for chaos to appear. 
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FIG. 22: The natural logarithms of the longest ellipsoid 
axis r vs. t for a strongly chaotic (20 + 1O)M0 configuration, 
for varying post-Newtonian terms. The initial spins (both 
maximal) are Si = (-0.935125, 0.329567, 0.130101) m? and 
S2 = (0.039523, -0.54303, -0.838783) mi while the initial 
radius is r = 4.318 Af, corresponding to = 240 Hz; 

the other initial conditions are fixed by the conditions for 
quasicircularity (Sec. UTTSl . The final time is 50 times the 
inspiral time, and the slopes of the lines are the Lyapunov 
exponents. In all cases, the Newtonian, IPN, 2PN, and spin- 
orbit terms are turned on, but some of the others may be 
turned off; from top to bottom (with colors, visible in elec- 
tronic versions of this chapter, noted parenthetically): all PN 
terms (black); 5iSi and S2S2 turned off (cyan); Si Si turned 
off (blue); 5*25*2 turned off (wiggly/red, nonchaotic); 5iSi, 
S2S2, and 5i52 turned off (straight/orange, nonchaotic). 



important result is that the spin-spin terms are crucial 
to the presence of chaos; when only Newtonian, IPN, 
2PN, 3PN and spin-orbit are turned on, the system is 
not chaotic. In the case illustrated in Fig. |^ the 5*15*2 
term by itself causes chaos: the presence of either the 
spin quadrupole term is irrelevant. This is not a general 
result: Fig. 1^ shows a case where the 5i5i quadrupole 
term apparently exerts a stabilizing influence: when only 
5252 is turned off, the system is nonchaotic, but if 5i5i 
is then turned off as well the system returns to chaotic 



behavior. 



C. Eccentric orbits 

Although quasicircular orbits represent the most likely 
source of gravitational waves from compact binaries de- 
tectable by ground-based observatories, some sources 
may consist of binaries with non neghgible eccentricity 
(such as produced, for example, by the Kozai mecha- 
nism ^J)- It is therefore potentially relevant to in- 
vestigate chaos for eccentric orbits. 



1. A chaotic eccentric orbit 

In agreement with 0, [S^, we find that the post- 
Newtonian equations can produce chaotic eccentric or- 
bits; one example appears in Fig. 1231 It must be em- 
phasized that, as in the case of quasicircular orbits, the 
pericenter of the orbit corresponding to Fig. [23 is at 
Tp = 5 M ; thus the two black holes are rather close to 
each other and higher order PN corrections, especially 
higher spin PN corrections, should also be taken into ac- 
count. Nevertheless, it is clear that the equations them- 
selves can produce chaotic solutions, even if those solu- 
tions have dubious physical relevance. 



2. Parameter variation 

In order to give a sense of the relative importance of 
various orbital and spin parameters for the presence of 
chaos, we take the system shown in Fig. and vary 
several parameters independently. In Fig. 1241 we show 
the effect on the dimensionless Lyapunov exponent A^insp 
of varying the pericenter. Since the various spin terms 
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FIG. 23: The natural logarithms of the longest ellipsoid 
axis r vs. t for a chaotic eccentric orbit. The slope of the 
hne gives Atinsp ~ 1.01, so that nearby trajectories diverge 
by a factor of e in approximately one inspiral time. The ini- 
tial conditions are Mi = 6, M2 = 4, X = (5,0,0), P = 
(0,0.61644,0.36160), Si = (-1.21570,-0.31859,0.81886), 
and S2 = (-0.15273,0.64525,-0.06902). For reference, we 
show the corresponding plot for a nonchaotic orbit (with the 
same initial conditions as above except for Si =0). 

(which make chaos possible) are decreasing functions of r, 
we might expect that the chaos is weaker or non-existent 
as the spin terms get smaller, and indeed this is the case: 
for the model system we consider here, there is no chaos 
for Tp > 6 M . Decreasing the eccentricity has a similar ef- 
fect (Fig.l25|l: more highly eccentric orbits are more likely 
to be chaotic, probably because the larger velocities lead 
to larger values of the nonlinear velocity-dependent terms 
in the Hamiltonian. Finally, varying the value of the spin 
parameter for one of the bodies (Fig. l26|l produces the ex- 
pected result: as the spin decreases, the chaos is generally 
suppressed or disappears altogether. From Figs. EIH^ 
it is clear that adding the third post-Newtonian term has 
only minor effects in general on the presence of chaos. 

3. A survey of eccentric orbits 

We undertake here a survey of eccentric orbits in an 
effort to understand the prevalence of chaos in these sys- 
tems. Unfortunately, in contrast to the quasicircular case 
(Sec. II V Bll , any survey of eccentric orbits is limited by 
the large number of parameters, which makes a compre- 
hensive survey impractical. Nevertheless, we have exam- 
ined thousands of initial conditions for a variety of masses 
and eccentricities, with special attention paid to systems 
that are realistic sources of gravitational radiation for 
ground-based detectors. 

We consider binaries with masses of (6 + ?>)Mq, (6 -I- 
4)M0, and (12-f 3)AfQ. Our choices for the eccentricities 
and pericenters are then guided by astrophysical consid- 
erations [27LI23 I: eccentricities are not larger than ~ 0.33, 
and we choose the orbital frequency at pericenter such 
that the corresponding (Newtonian) gravitational-wave 
frequencies lie in the frequency band of ground-based 
interferometers (for low eccentricities GW radiation is 
emitted mostly at one, two, and three times the orbital 



frequency). The specific values we consider here, which 
satisfy the conditions above, are e = 0.01,0.2,0.33 and 
/orb = 13.3,20,40,50,100 Hz. [Note that at eccentric- 
ity ~ 0.22 (~ 0.33) the amplitudes of the Newtonian 
gravity-wave signal for the first and third harmonics are 
17% (28%) and 40% (60%) relative to the second har- 
monic 30|.l The pericenter is then obtained from the 
formula 31] 



With the choices made for e and /orb, Eq- H4.ti|l takes 
values between 8M and 30 Af. Thus, the PN expansion 
is valid for the two compact bodies even at the point of 
closest approach. 

For each value of the binary masses and for each pair 
(e, /orb) we produce 500 randomly oriented maximal ini- 
tial spins and calculate the largest Lyapunov exponent as 
in Sec. IIVBI using a final integration time of ten times 
the inspiral time (as calculated from [l^). [For each ran- 
dom orientation we first fix the spin, and then use the 
method described in Sec. IIII Al to find the initial con- 
ditions; in particular, the PN Hamiltonian in Eq. I|3.3I) 
always includes the spin contributions.] In all cases con- 
sidered, we find no evidence of chaos — all the Lyapunov 
exponents are consistent with zero. 

V. CONCLUSIONS 

The dynamics of binary black holes, as modeled by 
the post-Newtonian equations, are significantly affected 
by the presence of spin. In particular, the addition of 
spin terms to the post-Newtonian equations of motion 
leads to significant changes in the orbital geometry and 
dynamical behavior of the solutions. The effects of the 
spin terms are particularly clear on quasicircular orbits, 
where the interaction terms quadratic in the spin cause 
deviations from perfectly spherical orbits. 

We find that, for quasicircular orbits, the presence of 
the interaction terms quadratic in the spins can lead to 
chaotic solutions, as indicated by positive Lyapunov ex- 
ponents. These exponents come in ±A pairs, a reflection 
of the Hamiltonian nature of the dynamics. We mea- 
sure the strength of the chaos by comparing the e-folding 
timescale for chaotic behavior (the inverse of the Lya- 
punov exponent) with the inspiral timescale. We find 
especially strong chaos for high-frequency/low-radius or- 
bits and high spins. However, in those cases the black 
holes are so close to each other that spin-spin induced 
oscillations in the radial separation are rather important, 
and the quasicircular initial conditions used in this pa- 
per may not correspond to widely separated quasicircu- 
lar orbits evolved adiabatically to low separations under 
gravitational radiation reaction. It would therefore be 
preferable to set initial conditions when the black holes 
are rather far apart (so that spin effects are negligible) 
and then evolve the system including radiation-reaction 
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FIG. 24: The dimensionless Lyapunov exponent Xtinsp as a function of pericenter for (a) all terms through 2PN and (b) all terms 
through 3PN. The initial spins are the same as those in Fig. 1231 with the other initial conditions fixed by the parameterization 
method described in Sec. lili Al Note that the Lyapunov exponent decreases with increasing pericenter as the spin coupling 
terms get smaller (consistent with the the results in 6;| ) . 
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FIG. 25: The dimensionless Lyapunov exponent Atinsp as a function of orbital eccentricity for (a) all terms through 2PN and 
(b) all terms through 3PN. Only eccentricities greater than around 0.6 show any evidence of chaos. The initial spins are the 
same as those in Fig. 1231 with the other initial conditions fixed by the parameterization method described in Sec. 1111 Al 



effects. The latter will be soon available in the Hamil- 
tonian framework, including spin couplings, and we plan 
to include them and investigate the presence of chaos in 
the near future. 

We build a survey of eccentric orbits which we believe 
is representative of binaries which can emit GWs in the 
LIGO/ VIRGO frequency band, have eccentricities justi- 
fied by astrophysical considerations, and whose dynamics 
can be safely described by the PN expansion. For this 
survey we do not find any chaos. We find chaotic solu- 
tions only for rather eccentric orbits (e ~ 0.9) with very 
low pericenters, which are not astrophysically motivated 
and for which higher order PN terms, especially higher- 
order spin couplings, should be consistently added in the 
ADM Hamiltonian. 

Since we find that chaotic behavior is due to spin- 
spin couplings and does not seem to be much affected 
by the non-spinning PN dynamics, we suspect that our 
results will not change qualitatively if the non-spinning 
ADM Hamiltonian were replaced by the Schwarzschild 
deformed effective-one-body Hamiltonian [s^ l33l| (which 
is a re-summation of the ADM Hamiltonian). By con- 
trast, we might expect differences if the Kerr-deformed 
effective-one-body Hamiltonian [TBj were used instead of 



the spinning ADM Hamiltonian. 

In conclusion, considering our surveys of quasicircular 
and eccentric orbits together, we find no chaos in any 
system for orbits that are of astrophysical interest for 
ground-based interferometers and which clearly satisfy 
the approximations required for the equations of motion 
to be physically valid at the post-Newtonian order con- 
sidered. 
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APPENDIX: NUMERICAL IMPLEMENTATION 

Our primary implementation of the post-Newtonian 
equations of motion is a collection of Mathematica pack- 
ages, relying on the native NDSolve function to effect 
numerical integrations. We implement the equations in 
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the parameterization method described in Sec. IIII Al The strength of the chaos decreases as the spin decreases, becoming 
indistinguishable from zero below around Si = 0.5 m^. 



a standard way to eliminate the mass variables, measur- 
ing lengths and times in terms of the total mass M and 
momenta in terms of the reduced mass fi. In these units, 
consistency then forces the angular momenta to be mea- 
sured in terms of /iM. 

We perform various tests to check the validity of 
the equations of motion and the numerical integrations. 
Most important, we check that the energy, angular mo- 
mentum, and spin magnitudes are all conserved by the 
integrations. These quantities are useful since they tend 
to be sensitive to mistakes in the equations of motion; we 
verify in all cases that the constants of the motion are 
conserved at a level consistent with the default accuracy 
goal (typically 10~^°), which gives us confidence that the 
equations are correct. We also check that orbital angular 
momentum and L • SoS are conserved through spin-orbit 
coupling, as discussed in Sec. Ill Bl Furthermore, we ver- 
ify that our implementation of the equations reproduces 
the Keplerian orbits of Newtonian gravitation (only New- 
tonian terms turned on), Lense-Thirring precession (ex- 
treme mass-ratio mi 3> m2 limit with Newtonian and 
SO terms turned on), and classical quadrupole preces- 
sion (Newtonian and S1S2 turned on). 

A second implementation of the equations of motion 
uses a Mathematica interface for a C-l — h numerical mr 

* The additional use of the freely available package Optimize .m |34| 
gives a factor of 5 increase in speed for the case at hand. 



tegrator. We use Mathematica to generate the required 
derivative expressions directly from the Hamiltonian, us- 
ing the built-in CForm function to convert to CH — h source 
code, so that the native Mathematica and C++ integra- 
tors automatically agree. ^ As in the case of the pure 
Mathematica integrator, we check the conservation of all 
the relevant quantities, in this case using a Bulirsch-Stoer 
integrator with an error goal of 10"^" (the variable eps 
from Numerical Recipes |35j). We also check several or- 
bits with a Runge-Kutta integrator with an error goal 
of 10"^; the resulting agreement with the Bulirsch-Stoer 
integrations verifies that our results are not specific to 
the choice of integration algorithm. 

As noted in Sec. lIV Al for the calculation of Lyapunov 
exponents we use the techniques and routines described 
in and especially Chapter 5 of 0|. In short, we use the 
deviation vector method for speed, but check select re- 
sults using the slower but more robust Jacobian method. 
We also use the Jacobian diagnostic described in 5] and 
[Tfij l to verify the correctness of the (rather complicated) 
Jacobian matrix (which in fact is generated by Math- 
ematica directly from the equations of motion). As a 
result, we are confident that our values for the Lyapunov 
exponents faithfully describe the true dynamics of the 
system. 
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